##################
##################
###DiD Analysis###
##################
##################
library(lfe)
library(sf)
library(spdep)
library(foreign)
library(doBy)
library(splm)
library(plm)
library(lmtest)
library(sandwich)
library(stargazer)
library(mgcv)

###Import data
#Set working library
setwd("~/OneDrive - Indiana University/FromGoogle/AfroGrid/Ghana/Code_6_4_25/data/")
#Read dataset
ghan.burk.dat <- read.csv("ghandatts.csv")


##############
###Main DiD###
##############
###Baseline
lm.1 <- felm(acled_mil_tot ~ lagacled_mil_tot+
               invbdist*y_2015_p|month+STATE_PROVINCE|0|gid,
             data=ghan.burk.dat)
summary(lm.1)

###Full
lm.2 <- felm(acled_mil_tot ~ lagacled_mil_tot+logNTL+logpop+logprec+t_c_avg+t+
               invbdist*y_2015_p|month+STATE_PROVINCE|0|gid,
             data=ghan.burk.dat)
summary(lm.2)

##Export to Latex
stargazer(lm.1, lm.2)

###Create marginal effects plot
##Call plot function
setwd("~/OneDrive - Indiana University/FromGoogle/AfroGrid/Ghana/Code_6_4_25/")
source("plots.R")
##Baseline model
setwd("~/OneDrive - Indiana University/FromGoogle/AfroGrid/Ghana/Code_6_4_25/plots/")
jpeg("figure2a.jpeg", width = 6, height = 6, units = 'in', res = 500)
felm_marginal_effects(regression_model = lm.1, 
                      data = ghan.burk.dat, 
                      treatment = "invbdist", moderator = "y_2015_p", 
                      treatment_translation = "inverse distance to border", 
                      moderator_translation = "Timing", 
                      dependent_variable_translation = "militia violence")
dev.off()

##Full model
setwd("~/OneDrive - Indiana University/FromGoogle/AfroGrid/Ghana/Code_6_4_25/plots/")
jpeg("figure2b.jpeg", width = 6, height = 6, units = 'in', res = 500)
felm_marginal_effects(regression_model = lm.2, 
                      data = ghan.burk.dat, 
                      treatment = "invbdist", moderator = "y_2015_p", 
                      treatment_translation = "inverse distance to border", 
                      moderator_translation = "Timing", 
                      dependent_variable_translation = "militia violence")
dev.off()


##########################
###Sensitivity analyses###
##########################

###Geospatial conflict models
#Balance the panel
afro.plm<- pdata.frame(ghan.burk.dat,index=c("gid","t"))
afro.plm.times<- make.pbalanced(afro.plm,balance.type = "shared.times")
afro.plm.balanced<- make.pbalanced(afro.plm.times,balance.type = "shared.individuals")
#Create gid indicators
gid<- data.frame(unique(afro.plm.balanced$gid))
gid$indiv<- 1:32
names(gid) <- c("gid","indiv")
##Read in the grid shape file
setwd("~/OneDrive - Indiana University/FromGoogle/AfroGrid/Ghana/Code_6_4_25/data/")
polygons.prio<- st_read("priogrid_cellshp/priogrid_cell.shp")
#Rearrange by gid
polygons.prio<- polygons.prio %>% dplyr::select(gid)
polygons.prio.slice<- subset(polygons.prio, gid %in% afro.plm.balanced$gid)
#Create grid indicator
rownames(polygons.prio.slice)<- 1:32
#Redesign the shapefiles as spatial tools
slice.sp<- as(polygons.prio.slice,"Spatial")
slice.nb<- poly2nb(slice.sp,queen=T)
# Row standardized
slice.listw<- nb2listw(slice.nb,style = "W",zero.policy = T)
#Redefine index variables as the first two columns
afro.plm.balanced<- merge(afro.plm.balanced,gid,by="gid")

####Sens. 1: Pooled standard errors by GID
spl.2.p <- spml(acled_mil_tot ~ lagacled_mil_tot+logNTL+logpop+logprec+t_c_avg+
                  invbdist*y_2015_p+factor(month)+factor(STATE_PROVINCE),
                data=afro.plm.balanced,listw=slice.listw,na.action=na.omit,
                index=c("indiv","t"),spatial.error = "b",
                lag=T,model="pooling")
summary(spl.2.p)

####Sens. 2: Random effects by GID
spl.2.r <- spml(acled_mil_tot ~ lagacled_mil_tot+logNTL+logpop+logprec+t_c_avg+
                  invbdist*y_2015_p+factor(month)+factor(STATE_PROVINCE),
                data=afro.plm.balanced,listw=slice.listw,na.action=na.omit,
                index=c("indiv","t"),spatial.error = "b",
                lag=T,model="random")
summary(spl.2.r)


######State/Province cluster SEs
####Sens. 3: Without STATE_PROVINCE FEs
lm.2.cy1 <- felm(acled_mil_tot ~ lagacled_mil_tot+logNTL+logpop+logprec+t_c_avg+t+
                   invbdist*y_2015_p|month|0|STATE_PROVINCE,
                 data=ghan.burk.dat)
summary(lm.2.cy1)

####Sens. 4: With STATE_PROVINCE FEs
lm.2.cy2 <- felm(acled_mil_tot ~ lagacled_mil_tot+logNTL+logpop+logprec+t_c_avg+t+
                   invbdist*y_2015_p|month+STATE_PROVINCE|0|STATE_PROVINCE,
                 data=ghan.burk.dat)
summary(lm.2.cy2)


####Sens. 5: Count models
#Create cubic time polynomials
ghan.burk.dat$t2 <- ghan.burk.dat$t^2
ghan.burk.dat$t3 <- ghan.burk.dat$t^3
#Rub model
glm.2 <- glm(acled_mil_tot ~ lagacled_mil_tot+logNTL+logpop+logprec+t_c_avg+t+t2+t3+
               invbdist*y_2015_p+factor(month)+factor(STATE_PROVINCE),
             data=ghan.burk.dat, family = "poisson")
#summary(glm.2)
glm.2.re <- coeftest(glm.2, vcov = vcovCL(glm.2, cluster = ~gid))

####Sens. 6: Extended environmental controls
lm.2.en <- felm(acled_mil_tot ~ lagacled_mil_tot+logNTL+logpop+logprec+t_c_avg+t+
                  NDVI_mean + spei + cropland_loose_portion+
                  invbdist*y_2015_p|month+STATE_PROVINCE|0|gid,
                data=ghan.burk.dat)
summary(lm.2.en)


####Sens. 7: Other conflict controls 
lm.2.co <- felm(acled_mil_tot ~ lagacled_mil_tot+logNTL+logpop+logprec+t_c_avg+t+
                  acled_battle_state + acled_battle_rebel + 
                  invbdist*y_2015_p|month+STATE_PROVINCE|0|gid,
                data=ghan.burk.dat)
summary(lm.2.co)

####Sens. 8: Extended environmental + conflict controls
lm.2.coen <- felm(acled_mil_tot ~ lagacled_mil_tot+logNTL+logpop+logprec+t_c_avg+t+
                    acled_battle_state + acled_battle_rebel + 
                    NDVI_mean + spei + cropland_loose_portion+
                    invbdist*y_2015_p|month+STATE_PROVINCE|0|gid,
                  data=ghan.burk.dat)
summary(lm.2.coen)

####Sens. 9: Alternative temporal controls (comparable to SCs)
lm.2.tc <- lm(acled_mil_tot ~ lagacled_mil_tot+logNTL+logpop+logprec+t_c_avg+t+
                invbdist*y_2015_p+factor(month)+t*factor(STATE_PROVINCE),
              data=ghan.burk.dat)
#summary(lm.2.tc)
lm.tc <- coeftest(lm.2.tc, vcov = vcovCL(lm.2.tc, cluster = ~gid))


####Sens. 10: GMM
#Convert data to pgmm formal
m.pgmm <- pdata.frame(ghan.burk.dat, index=c("gid", "t"))

#set seed
set.seed(50)

##Run model
gmm.2 <- pgmm(acled_mil_tot ~ lagacled_mil_tot+logNTL+logpop+logprec+t_c_avg+
                invbdist*y_2015_p| 
                lag(acled_mil_tot, 2:5),
              data = m.pgmm, effect = "individual", transformation= "ld", model = "onestep")
summary(gmm.2)


####Sens. 11: Time-varying unobserved confounders (Province month)
#create a quantitative version of STATE_PROVINCE
ghan.burk.dat$prov_num <- ifelse(ghan.burk.dat$STATE_PROVINCE=="Northern",1,
                                 ifelse(ghan.burk.dat$STATE_PROVINCE=="Upper West",2,3))
ghan.burk.dat$STATE_PROVINCE_month <- ghan.burk.dat$prov_num*ghan.burk.dat$month
#Run model
lm.2.pm <- felm(acled_mil_tot ~ lagacled_mil_tot+logNTL+logpop+logprec+t_c_avg+t+
                  invbdist*y_2015_p|STATE_PROVINCE_month|0|gid,
                data=ghan.burk.dat)
summary(lm.2.pm)

####Sens. 12: Time-varying unobserved confounders (Province month)
#Create interaction
ghan.burk.dat$STATE_PROVINCE_year <- ghan.burk.dat$prov_num*ghan.burk.dat$year
#run model
lm.2.py <- felm(acled_mil_tot ~ lagacled_mil_tot+logNTL+logpop+logprec+t_c_avg+t+
                  invbdist*y_2015_p|month+STATE_PROVINCE_year|0|gid,
                data=ghan.burk.dat)
summary(lm.2.py)

####Sens. 13: Nonlinear dynamics (Additive Model)
gam2<- gam(acled_mil_tot ~ s(invbdist, by = y_2015_p)+
                   lagacled_mil_tot+logNTL+logpop+logprec+t_c_avg+t+
                   factor(STATE_PROVINCE) + factor(month), data=ghan.burk.dat)
summary(gam2)
#Plot the smoothed trend
setwd("~/OneDrive - Indiana University/FromGoogle/AfroGrid/Ghana/Code_6_4_25/plots/")
jpeg("figureA2.jpeg", width = 6, height = 6, units = 'in', res = 500)
plot(gam2)
dev.off()
#summary(1/ghan.burk.dat$invbdist) shows it increases and peaks at around 70KM from the border


####Export to Latex
#Table A1
stargazer( lm.2.cy1, lm.2.cy2, glm.2.re)
#Table A2
stargazer(lm.2.en, lm.2.co, lm.2.coen, lm.tc, gmm.2)
stargazer(lm.2.en, lm.2.co, lm.2.coen, lm.2.tc, gmm.2)
#Table A3
stargazer( lm.2.pm, lm.2.py, gam2)
